Systematic Semiclassical Expansion 
for Harmonically Trapped Ideal Bose Gases 
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Using a field-theoretic approach, we systematically generalize the usual semiclassical approxi- 
mation for a harmonically trapped ideal Bose gas in such a way that its range of applicability 
is essentially extended. With this we can analytically calculate thermodynamic properties even 
for small particle numbers. In particular, it now becomes possible to determine the critical tem- 
perature as well as the temperature dependence of both heat capacity and condensate fraction in 
low-dimensional traps, where the standard semiclassical approximation is not even applicable. 
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I. INTRODUCTION 

The field of ultracold Bose gases attains at present a lot 
of attention due to an improved experimental accessibil- 
ity within the last decade. Many different theoretical ap- 
proaches are used to treat these trapped dilute quantum 
gases. Although isolated Bose gases should, in principle, 
be described within the micro-canonical ensemble, one 
commonly applies the technically more efficient canonical 
or grand-canonical descriptions [1-20]. This is justified 
as experiments often use a large number of bosons. It is a 
common belief in quantum statistics that, at least in the 
thermodynamic limit N— >oo, all ensembles should con- 
verge to one and the same result. However, we note that 
some peculiar exceptions are known for particle count- 
ing statistics as discussed, for instance, in Refs. [20-22]. 
From a theoretical point of view, the grand-canonical en- 
semble has the advantage that it provides an analytical 
description, whereas the canonical approach is limited to 
numerical results for moderate particle numbers. As ex- 
periments with ultracold Bose gases are always realised 
with a finite number of particles, the fundamental ques- 
tion arises how to study finite-size effects for the ther- 
modynamic properties of trapped Bose gases most effi- 
ciently. 

To analyse this problem systematically, we introduce 
and compare two different approaches. In Sect. II we 
briefly rederive the well-known grand-canonical descrip- 
tion without using the order parameter concept for a har- 
monically confined ideal Bose gas. The thermodynamic 
properties can only be calculated numerically in this the- 
ory. In Sect. Ill we introduce another grand-canonical 
description of the trapped Bose gas which is analytical 
as it relies on the order parameter concept. Introduc- 
ing an order parameter is an essential approximation for 
finite systems and leads to different results for the re- 
spective thermodynamic quantities compared to the the- 
ory without order parameter. However these differences 
vanish in the thermodynamic limit and turn out to be 
negligibly small for experimentally realistic system sizes. 



Moreover, generalizing a formalism developed in Ref. [23, 
Appendix 7A] , this approach extends the usual semiclas- 
sical approximation [1, 2, 13, 16] to a systematic semi- 
classical expansion which yields yet unknown analytical 
results for the thermodynamic quantities in the super- 
fluid phase. In particular, we will calculate the critical 
temperature T c as well as the temperature dependence of 
the condensate fraction No/N and the heat capacity Cy 
for D=l, 2, 3 dimensions up to the order of the semiclas- 
sical expansion which was not accessible before. 



II. APPROACH WITHOUT ORDER 
PARAMETER 



We start with briefly rederiving the well-known grand- 
canonical description of an ideal Bose gas. The general 
expression for the grand-canonical potential of an ideal 
Bose gas is given by [24] 
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where n describes the one-particle quantum numbers. 
Here E n , /3=l/fceT, and fi denote the energy levels of the 
system, the inverse temperature, and the chemical poten- 
tial, respectively. We specify Eq. (1) for the case of an 
ideal Bose gas which is trapped in an isotropic harmonic 
potential of the form I/(x)=Mw 2 x 2 /2, where M and ui 
denote the mass of a bosonic particle and the trap fre- 
quency, respectively. The one-particle energy eigenvalues 
ofaharmonic oscillator are E n — fiw(ni + . . .+rirj + D/2) 
in D dimensions, so the grand-canonical potential (1) 
specifies to 
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For numerical calculations it turns out to be useful to 
follow Ref. [25] and reexpress (2) by using the polyloga- 
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FIG. 1: Dimensionless chemical potential (fi—Eo)/hu versus 
T/T c in D—3 dimensions determined from (6) for TV=10 3 par- 
ticles (dashed line) and JV=10 4 (solid line) for lj—2tx-4Q Hz. 
The critical temperatures are T c =16.75 nK for TV=10 and 
T c =37.56 nK for TV=10 4 , respectively. 
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FIG. 2: Heat capacity (8) versus T/T c in D=3 dimensions for 
TV=10 4 (dashed line) and TV=10 5 (solid line) for a;=27r-40 Hz. 
The critical temperatures are T c =37.56 nK for TV=10 4 and 
T c =82.46 nK for TV=10 5 , respectively. 
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and the identity 
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With N=—dtt/dfj, one gets for the particle number equa- 
tion 

N = E ( m+ ^ _1 ) Co (e^"-^-" 1 ^) , (6) 

which can be solved numerically for the chemical poten- 
tial /i once the particle number TV is given. In this ap- 
proach we consider the resulting condensate fraction 
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internal energy U=fl+TS+^iN and yields 
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One can see in Fig. 2 that Cy in D~'i dimensions has 
its maximum at TkT c and vanishes exponentially fast 
in the limit T|0. The heat capacity Cy gets larger for 
T<T C and smaller for T>T C if the particle number TV 
is increased. Moreover, one obtains that the slope at 
T=T C increases as well for larger TV, but does not di- 
verge. With this we have shown exemplarily that phase 
transitions do not occur in finite systems. However, \x 
and Cy for fixed T/T c seem to tend towards a limit for 
large particle numbers TV and thus, we expect the emer- 
gence of a sharp phase transition in the thermodynamic 
limit TV^oo. In the next section we use this observation 
as a motivation to introduce an analytical approach for 
describing a trapped Bose gas with the help of an order 
parameter. 



III. APPROACH WITH ORDER PARAMETER 



as a function of the temperature and define the point 
where the curvature is maximal as the critical tempera- 
ture T c . Figure 1 shows the chemical potential fi in D=3 
dimensions determined from (6) as a function of temper- 
ature T for a given particle number TV. One observes 
that fx remains smaller than the ground-state energy Eo 
for all temperatures and approaches E$ in the limit TJ.0. 
Moreover, one can see that (fi— E ) / fku gets smaller for 
T<T C if the particle number TV is increased. 
The heat capacity Cy=dU/dT\v,N is derived from the 



We start with the functional integral approach to 
the grand-canonical partition function of a harmonically 
trapped ideal Bose gas [26] 

Z=<j>V^*ViPe- A[r -' 4 ' ]/h , (9) 

where one integrates over all possible bosonic Schrodingcr 
fields •?/?* (x, t), ^(^i 7 ") which are periodic in imaginary 
time r with period h/3. The Euclidean action A[ip*,ip] 
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reads 
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We evaluate the functional integral by using the back- 
ground method [27, 28]. To this end we divide the fields 
iP*(x,t), V(x, t) into field expectation values 5'*(x, r), 
^(x, t), which we identify later on with the macro- 
scopic occupation of the ground state, and fluctuations 
t), ^(x,t): 



V>*(x,t) = **(x,t) +#*(x,t) 
V>(x,t) = *(x,r) +5V(x,r). 



(11) 



Note that field expectation values and fluctuations have 
to satisfy the condition [29] 



d D x <J-*(x,t)£V(x,t) =0. 
Using (11) together with (9) we arrive at 
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Z = e- A ^-^' h j>V5rV5^e- A ^-^l h . (13) 

Now we decompose the fluctuations Sip(x, r) into the one- 
particle eigenstates Vn( x ) of the system and apply an 
additional Matsubara decomposition: 
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with the Matsubara frequencies Lo m =2irm/h(3. Note that 
we explicitly do not sum over the ground state of the 
system in (14), as we have to satisfy condition (12). With 
this, the measure of the functional integration (13) turns 
into 
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The integration over the expansion coefficients c* m , c n ^ m 
in (15) is now Gaussian and can be performed. The ef- 
fective action is then found by applying the logarithm to 
the partition function: T[**, *]= - (log Z) /f3. With this 
one gets 
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This effective action yields the grand-canonical potential 
fl if it is evaluated for extremised field expectation values: 



fJ=r[ , !l>*, <3> c ]. An extremization of (16) with respect to 
**(x, t) leads to 
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This equation has the eigenstates V"n(x) with 
/i=£' n as non-trivial solutions which are periodic 
in imaginary time. We choose the ground state 
ip (x)=(Muj/irh) D/4 cxp (-Mwx 2 /2fi) to be the physi- 
cally meaningful solution and normalise 5 , c (x, r) to the 
number of atoms in the ground state No'. 

* c (x,r)-v / ^o^o(x). (18) 
As (17) and (18) lead to the algebraic equation 
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(19) 



we obtain two different phases. In the gas phase we have 
A =0, whereas in the superfluid phase with A o ^0 the 
chemical potential /x must be equal to the ground-state 
energy E . The critical temperature T c occurs at the 
borderline between both phases, so it follows from the 
particle number equation by setting both N = and 
/i = Eq. To this end we combine (16) and (18) and get 
for the grand-canonical potential 
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The first term represents the contribution of the macro- 
scopically occupied ground state of the system, whereas 
the second term describes the thermal contributions of all 
excited states. This should be compared with (2) where 
the ground state is treated like all other states. 



IV. THERMODYNAMICS 

It is now possible to derive all thermodynamic quanti- 
ties from (20) within this framework. The particle num- 
ber equation N——dfl/dfi reads 



N = N a +J2< 
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The internal energy follows from the Legendre transfor- 
mation U=n+TS+nN: 



U = E a N a +Dhuj V e^- £ °) fc { 1 
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Both the particle number N and the internal energy U 
can be expressed by the auxiliary functions 
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The particle number (21) turns into 

N = N + (f ) W - »'> b > D ) > ( 24 ) 

where we have used the abbreviations b~hw/3 and 
fi'=((i—E )/huj. For the internal energy (22) we get cor- 
respondingly 



U = E N + Dhuj 



+I(l-li',b,D+l) 



(25) 



With this it is possible to calculate useful analytic ap- 
proximations for the critical temperature T c , as well as 
the temperature dependence of both the condensate frac- 
tion No/N and the heat capacity Cy for different num- 
bers of spatial dimensions D. To this end we assume 
that the difference ftuj between different energy levels of 
the harmonic oscillator is small compared to the average 
thermal energy k B T . This yields the semiclassical condi- 
tion 0<6=fiw/3<§Cl, which is well fulfilled for present-day 
experiments. 

In order to describe the supcrfluid phase we apply 
the limit ^Eq to the particle number equation (24) for 
D=l,2,3 and use the formulas (A13)-(A15) which are 
derived in the Appendix within the semiclassical approx- 
imation. Furthermore, we expand the polylogarithms for 
small b—hjj(3 by using the Robinson formula [30] 
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Thus, we get for T<T C 
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Here (,{z) denotes the Riemann zeta function and 
7=0.5772. . . is Euler's constant. Surprisingly, our results 
(27)-(29) coincide with the findings of two different ap- 
proaches which are reviewed and compared in Rcf. [12]. 
The first one is a master equation approach to canonical 
condensate statistics which is based on an analogy to the 
laser phase transtion [10, 19]. It yields accurate results 
even for small systems, is valid for all temperature, but 
is partly numerical. The second one, which is entirely 
based on considering the particle-number distribution, is 
fully analytical, but is limited to temperatures below T c 
[11, 15]. ' 
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FIG. 3: Finite-size corrections (T c -T c (0) )/T c (0) taken from 
Eq. (32) versus iV~ 1/3 for D=3 dimensions. The dotted 
(dashed) line includes the first (the first two) finite-size cor- 
rection^). The bullets correspond to the critical temperature 
which is numerically determined from (6) and (7) of the the- 
ory without order parameter. 



A. Critical Temperature 

Setting Ao=0 in (27)-(29), one obtains for the critical 
temperature 
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(32) 



In D=\ dimension the critical temperature T c follows 
from numerically solving the implicit Eq. (30). Higher 
corrections to (30) are of the order 0(N^ 1 ), so they 
are small and can be neglected. Note that our result 
(30) differs slightly from the corresponding finding of 
Refs. [2, 18]. For D=2,3 one defines the thermodynamic 
limit N^>oo in such a way that the leading order of (31) 

and (32), i. e. T c (0) ^hwN 1 / /k B C(3) 1/D , remains con- 
stant. Thus, when the particle number N is sent to in- 
finity the trap frequency w has to approach zero in such 
a way that the product Nuj d is kept fixed. In this case 
the quantity —huj (3^ is, indeed, small as has been 
assumed above. Higher orders of (31) and (32) are called 
finite-size corrections. We note that additional nontrivial 
logarithmic dependences on the particle number N occur 
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which do not follow from the standard semiclassical ap- 
proximation [1, 2, 13, 16]. In Fig. 3 the critical temper- 
ature in D=3 dimensions is plotted up to the first and 
the second order for different particle numbers N and 
compared with the corresponding finite-size corrections 
which are obtained numerically from the theory with- 
out order parameter from (6) and (7). Combining the 
first and second finite-size corrections from (32) yields 
a better agreement with the theory without order pa- 
rameter than the first finite-size correction alone. Fur- 
thermore, we read off from Fig. 3 that, for particle num- 
bers larger than about 7V=10 3 , the analytic formula (32) 
yields values for the critical temperature which coincide 
with the corresponding results of the theory without or- 
der parameter for all practical purposes. Thus, although 
introducing an order parameter for studying finite-size 
effects represents an essential approximation, its findings 
do not differ from the results of the theory without an 
order parameter for experimentally realistic system sizes. 
Note that the first correction of (32) was already found 
some time ago [1, 2], whereas the second correction has 
only recently been found [23, Appendix 7A]. 
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FIG. 4: Condensate fraction No/N versus temperature T for 
a) D—l dimension and N=W 4 and for b) D=2 dimensions 
and N=10 3 for lu=2tt-40 Hz. The solid lines describe Eq. (7) 
from the theory without order parameter, the dashed lines 
correspond to the leading orders of (33), (34) and the dotted 
line also includes the first finite-size correction. 



already for 7V=10 3 particles as shown in Fig. 4 b). Note 
that Eq. (34) including the first finite-size correction 
yields a better agreement with the theory without or- 
der parameter than the already known leading order of 
(34). 
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C. Heat Capacity 



The new feature of our approach is that it is applicable 
in the whole temperature regime. For instance, one gets 
for the condensate fraction N /N from (27)-(32) 
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The heat capacity Cy can be found from (25) by us- 
ing the relation Cy—dUjdT\ y N . One has to take into 
account two different regimes. For T<T C the chemical po- 
tential [i is fixed and the number of atoms in the ground 
state N depends on temperature. With this we get 
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Figure 4 shows how the approximation to describe a fi- 
nite system with an order parameter deviates from the 
original grand-canonical approach without order param- 
eter. In Fig. 4 a) one observes in D—l dimension that 
the approximation introduced by the decomposition (11) 
is still noticeable at TkT c for 7V=10 4 particles. For D=2 
dimensions, however, this effect becomes much smaller 



For T>T Cl on the other hand, Nq vanishes and \i depends 
on temperature. Thus, the heat capacity for T>T C still 
depends explicitly on T and //. However, we can analyt- 
ically work out the limit TJ,T C and obtain 
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Thus, the heat capacity has a discontinuity at T=T C in 
£>=2,3 which follows from (37), (38), (40), and (41): 
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In 73=3 dimensions the discontinuity ACy remains fi- 
nite in the thermodynamic limit TV— >oo and gets smaller 
for finite systems. This can be seen in Fig. 5, where 
the heat capacity taken from (38) and (41) including the 
first finite-size correction are compared with the thermo- 
dynamic limit. In D=2 dimensions one has a disconti- 
nuity for finite systems which vanishes in the thermody- 
namic limit TV— »oo. Finally, in D=l dimension there is 
no discontinuity ACy at T=T C at all in leading order in 
agreement with the findings of Rcf. [18]. Note that the 
leading contribution in the heat capacity discontinuities 
ACy differs slightly from the findings in Rcf. [18]. 

Note that Eqs. (36)-(38) also provide systematic semi- 
classical expansions for both the internal energy U and 
the entropy S by using the thermodynamic relations 
C v =dU/dT\ ViN =TdS/dT\ ViN [24]. However, one does 
not need to work out the limit T^T C separately above 
and below the critical point as both the internal energy 
U and the entropy S are continuous at T=T C . 




T[nK] 



FIG. 5: Heat capacity Cv/Nks versus Temperature T in 
D=3 dimensions for w=2tt-40 Hz and iV=10 4 . The dotted 
(dashed) line is taken from Eqs. (38), (41) and includes the 
leading order (the first two leading orders). The solid line 
corresponds to (8) from the theory without order parameter. 



V. CONCLUSIONS 

In this paper we have extended the usual scmiclas- 
sical expansion [1, 2, 13, 16] for harmonically confined 
ideal Bose gases. With this we have derived orders of 
the semiclassical expansion which have not yet been ac- 
cessible using standard semiclassical approaches. This 
has been shown exemplarily for the critical tempera- 
ture T c as well as the temperature dependence of the 
condensate fraction Nq/N and the heat capacity Cy in 
D=l,2, 3 dimensions. It would be straight-forward to 
generalize our findings to anisotropic harmonic trapping 
potentials, which are used in many experiments to study 
Bose-Einstein condensation in D = 2, 3. Finally, we note 
that it would be quite instructive to clarify the connection 
between our field-theoretic approach towards a system- 
atic semiclassical expansion and the statistical approach 
of Refs. [9, 11, 12, 14, 15] as both seem to be related (see 
the remark below Eq. (29)). However, this relation is by 
no means obvious as both the grand-canonical and the 
canonical approach rely on certain approximations. 
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APPENDIX A: SEMICLASSICAL 
APPROXIMATION 

Some thermodynamic properties of the ideal Bose gas 
in a harmonic trap are expressable in terms of the series 
(23). Generalizing an approach of Ref. [23, Appendix 7A] 
we work out a systematic semiclassical approximation of 
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(23) which is valid for 0<6<1. 

1. Euler-MacLaurin Formula 

In the limit of small b it is suggestive to approximate 
(23) by the Euler-MacLaurin formula 



fe=0 



dk f(k) + 



f(0) + f{g) 



+ ... . 
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However, this is not directly possible, as the integral 
would diverge if the series starts at k=0. One way to 
avoid this divergency is to subtract all divergent terms 
in (23) before replacing the series by an integral. To this 
end we expand the denominator of (23) for small b 



D-l 



(1 



\D E 



Ci(D) 



- bk) D {bk) D- 



J+O(b ) 



(A2) 



where C Q (D)=1, Ci(D)=D/2, C 2 {D)=D(3D-1)/2A, 
Cs(D)=D 2 (D— 1)/48, ... are the respective expansion 
coefficients. Afterwards we subtract this expansion from 
(23) and add it again: 



I(A,b,D) = J2 e 



bAk 



k=0 



Ci(D) 



oo D-l 

EE 

fc=i i=o 



(1-e- 
Ci(D) e 



-bk\D 



D-l 

S \bk) D - 
1=0 y ' 



-bAk 



(bk) 



D-l 



+ O{b ). (A3) 



As we have extended the summation in the first line from 
k=l to fc=0, we have obtained an error of the order O(b ). 
Finally we apply (Al) to the first line in (A3), where 
higher corrections can be neglected as they are as well of 
the order O(b ). With this we arrive at 



I(A,b,D) 



dk e 



-bAk 



D-l 



cm 

(l-e- bk ) D £jj (bk) D - 1 



d-E 



oo D 



2.2. 

k=l 1=0 



I Ci(D) er bAk 



■ O(b ) . (A4) 



{bk) D ~ l 

Both terms in (A4) are treated as follows: 

• We substitute x=e~ bk in the first integral of (A4). 
One immediately observes that the contribution of 
this integral is of the order 0{b~ r ). Moreover one 
can reexpress the first integral in terms of Gamma 
functions by using the Beta function [31, (8.380)] 



B(x,y) = f dtf-^l-t)"- 1 = 
Jo 



T(x + y) 



(A5) 



The remaining integrals in (A4) can directly be 
evaluated using the integral representation of the 
Gamma function [31, (8.310)]: 



f 

Jo 



dkk 



x-l 



e - ka = 



T(x) 



(A6) 



• The series in the second line of (A4) can be ex- 
pressed in terms of the polylogarithms (3). 

With this we obtain eventually 



I{A,b,D)=\ 



D-l 



]T Q(D) T(l+l-D) 



xA 



D-l-l 



T(1-D)T(A) 
T(l+A-D) ' 

E^CD-«(e-^)+0(6P). (A7) 



1=0 



2. Dimensional Regularisation 

As the Gamma functions in (A7) are divergent for in- 
teger dimension D, we apply dimensional regularization 
[32, 33]. To this end we set D=d-e with d=l, 2, 3, 4 . . . 
and consider the limit e^O. This leads to 



I(A,b,d) 



lim — 

e^O b 



T(l-d+e)T(A) 



T(l+A-d + e) 
^2Ci(d-e) T(l+l-d+e)A 



d-l 



d-l-l-. 



1=0 

d-l 



1=0 



+O{b ) . (A8) 



In order to evaluate the limit e^O we use the following 
expansion for the Gamma function 



r(-n+e) 



(-1)' 



-+M^+l)+0{e 1 ) 



(A9) 



with n=l,2,3, ... and e>0 [32, (8D.24)]. The digamma 
function ipo(z) is defined as ^po{z)~V '(z) /T(z). 
It satisfies the recursion formula [32, (8D.6)] 
V>o(2)=l/(2-l)+^o(2-l), where 7 =-^ (l)=0.5772 . . . 
denotes Euler's constant. Moreover, we use the 
approximation x £ —l+eliix+ . . . and expand the co- 
efficients Ci{d-e)^cl a \d)-Cf } {d)e+0{e 2 ) for small 
e, where cf ] (d)=C,(d), c£\d)=0, C{ 1) (d)=l/2, 

cf } (d)=(6d-l)/24, c£\d)=d(3d-2)/48, ... are the 
respective expansion coefficients. With this (A8) reduces 
to 



I(A, b, d)= lim ^ — I 

£^o b 



1 



(d-l)' 



d-l 



i=i 



— ipo^+A-d) +ip (d) 



E 

(=0 

+ E (^W )- 

(=0 



(-im (0) (^) 
(d-i-iy. 



A 



d-l-l 



(A10) 



Equation (A10) can be further simplified if one uses the 
identity 



(iijin^-o-E 

v ; 1=1 1=0 



1 (-lyc^d)^-,.! 



(d-i-iy. 



(AH) 



/(A ) 6,2) = ic2(e- A6 )+^-l^-l)[l» l-''o(.l-l)] 



+^+Ci( e - Ab )}+o(&°), 



(A14) 



which can be proven by complete induction. Thus, the 
terms in (A10), which are proportional to 1/e, cancel 
and we obtain a finite result in the limit e^O. Using the 
recursion formula of the digamma function one arrives at 



1=0 



{d-i-iy. 

cl°\d) 



d-l 



+ E 1 



m=d—l 



d-l n (0) 



+ 2^ w- 1 



-A6\ 



+0(6°) . 



(A12) 



Assuming 0<6<§;1 we read off that Eq. (A12) is, indeed, 
a good approximation for (23) as the error is of the or- 
der O(b ). For different dimensions d=l,2,3,4 one gets 
explicitly 

/(A, 6,1) = 1 [ln^-Vo(^)+Ci(e- A6 )] +O(ft ) ,(A13) 



'( ^'•■') = ^C3(e- A6 )+^C2(e- A6 )+ 7 ^-(,l-l) 

5 
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x (A-2) [In A-^(A-2)] --A + £ + Ci( e - Ab ) 



+0(6°) , 



(A15) 



^.M) = ^Ue- Ab )+y 3 Cs(e- Ab )+^- Ab ) 
x(A-3) [In A- Vo (^4-3)] |+O(6 ) . (A16) 



Note that the respective polylogarithmic functions in 
(A13)-(A16) have to be evaluated for small b by using 
the Robinson formula (26). 
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